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A thin thread of viscous fluid falhng onto a moving belt generates a surprising variety 
of patterns depending on the belt speed, fall height, flow rate, and fluid properties. 
Here we simulate this experiment numerically using the Discrete Viscous Threads 
method that can predict the non-steady dynamics of thin viscous filaments, captur- 
ing the combined effects of inertia, stretching, bending and twisting. Our simulations 
successfully reproduce nine out of ten different patterns previously seen in the labo- 
ratory, and agree closely with the experimental phase diagram of Morris et al. (2008). 
We propose a new classification of the patterns based on the Fourier spectra of the 
longitudinal and transverse motion of the point of contact of the thread with the belt. 
These frequencies appear to be locked in most cases to simple ratios of the frequency 
Qc of steady coiling obtained in the limit of zero belt speed. In particular the in- 
triguing 'alternating loops' pattern is produced by combining the first five multiples 
of fie/3. 
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I. INSTABILITIES OF VISCOUS THREADS 
A. Introduction 

A thin stream or jet of liquid falling onto a fixed surface is one of the simplest situations 
in fiuid mechanics, yet it can generate a remarkable range of phenomena. Fast jets produce 
hydraulic jumps, which are circular when the viscosity is very Ioa^P^ and polygonal when 
it is somewhat higher^. Thin streams of very viscous fiuid can exhibit steady coiling or 
rotatory folding, and under some conditions coiling produces spiral waves of air bubbles in 
the thin fiuid layer spreading over the surfac^. Finally, thin streams of non- Newtonian fiuid 
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FIG. 1. Configuration of the fluid meclianical sewing macliine. Newtonian fluid with constant 
density viscosity z/, and surface tension coefficient 7 is ejected at a volumetric rate Q through 
a nozzle of diameter d at a height H above a belt moving in its own plane at constant speed V . 
The position of the contact point between the thread and the belt is x(t). Lateral advection of the 
contact point motion at speed V creates complex 'stitch' patterns on the belt. In inset: geometry 
of an element of the thread, showing the orthonormal triad of basis vectors di, d2, and ds as a 
function of the Lagrangian coordinate S along the center-line. 

can exhibit the Kaye ("leaping shampoo") effect in which the stream rebounds episodically 
from the pile of previously deposited fiuicP. 

A further degree of complexity is introduced if the source of the jet and the surface onto 
which it falls are in relative motion. This is the case when a home cook lays down 'squiggles' 
of icing or molten chocolate on a cake, or when an artist lets paint dribble onto a canvas from 
a moving paintbrush, a technique used to great effect by Jackson PoUockP. An analogous 
situation involving many interacting jets is the 'spunbonding' process of non-woven fabric 
production, in which thousands of threads of molten polymer solidify and become entangled 
as they fall onto a moving belt, producing a fabric with a random texture. 

Here we use a numerical approach to study an idealized model for these processes: the 
continuous fall of a single thread of viscous fluid onto a belt moving with a constant velocity 
in its own plane (Figure [T]). This system was flrst studied experimentally by Chiu- Webster 
and Listei'^ (henceforth CWL), who called it the 'fluid mechanical sewing machine' on ac- 
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count of the stitch-like patterns traced on the belt by the thread. The complexity and 
diversity of these patterns testifies to a rich nonlinear dynamics and bifurcation structure. 
The appeal of the system is further increased by the theoretical and numerical challenges 
involved in modeling it. 

In CWL's experiments, viscous fluid (corn syrup) with density p, surface tension coeffi- 
cient 7 and viscosity v was ejected at a volumetric rate Q from a vertical nozzle of diameter 
d, from which it fell a distance H onto a belt moving at speed V (fig. [T|. The experiments 
were conducted by varying V and H for several different combinations of values of (i, Q and 
V. When V greatly exceeded a fall height-dependent critical value Vb(i^), the fluid thread 
had the form of a steady dragged catenary. As the belt speed was decreased towards VJ,, the 
lowermost part of the thread evolved into a backward-facing 'heel', which became unstable 
to periodic meandering when V = Vb- Further decrease of the belt speed led to a series of 
bifurcations to more complex patterns (flg. [sj, ending with the establishment of steady coil- 
ing for V = 0. CWL successfully predicted the shape of the steady dragged catenary using a 
'viscous string' model that neglected bending stresses. However, this solution ceases to exist 
when the extensional axial stress at the bottom of the thread becomes zero, corresponding 
to the incipient formation of a heel in which the axial stress is compressional. Because a 
state of axial compressive stress is a necessary condition for the buckling of a slender bod}^^, 
CWL interpreted the onset of meandering as a buckling instability of the heel. 

Ribe et al.^^ carried out a numerical linear stability analysis of the dragged catenary state 
to predict the critical belt speed Vb and the frequency Ub for the onset of meandering, using 
a more complete 'viscous rod' theory incorporating bending and twisting of the fllament. 
The numerical predictions of Vb and Ub thereby obtained agree closely with the experimental 
measurements oM Ribe et al.'^ also documented a close correspondence between incipient 
meandering and finite-amplitude steady coiling on a motionless {V = 0) surface, such that Ub 
is nearly identical to the steady coiling frequency Qc for any given fall height H. Moreover, 
the critical belt speed Vb{H) is nearly identical to the vertical (free-fall) speed f// of the fiuid 
at the bottom of the thread, indicating that meandering sets in when the belt is no longer 
moving fast enough to carry away in a straight line the fiuid falling onto it. 

More extensive experiments were conducted by Morris et al.^^, using a carefully engi- 
neered apparatus with silicone oil as the working fiuid for better stability and reproducibility. 
They determined a complete phase diagram for the patterns as a function of H and V for 
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a particular set of values of the fluid viscosity v and the flow rate Q. They showed that 
the observed amplitude of meandering as a function of the belt speed is consistent with a 
Hopf bifurcation, and proposed a simple model to predict it based on the hypothesis that 
the contact point moves at constant speed relative to the belt. Finally, they proposed a 
generic set of amplitude equations which they used to characterize the alternating loops 
(which they called 'figure-of-eight') and translated coiling patterns. 

Most recently, Blount and Listei'^ performed a detailed asymptotic analysis of a slender 
dragged viscous thread, focussing on the structure of the heel. They showed that the lower- 
most part of the thread can exhibit three distinct dynamical regimes depending on whether 
the belt speed is greater than, nearly equal to, or less than the free-fall speed f/j. Their 
asymptotic stability analysis of these steady states indicates that meandering sets in when 
the horizontal reaction force at the belt begins to be slightly against the direction of belt 
motion, corresponding to the heel 'losing its balance'. 

As the above summary indicates, our current theoretical understanding of the fluid- 
mechanical sewing machine is essentially limited to the initial bifurcation from the steady 
dragged configuration to meandering. In this paper, we push forward into the fully nonlin- 
ear regime with the help of a new computational algorithmP'^ that permits for the first 
time robust numerical modeling of arbitrary non-stationary dynamics of viscous threads. 
After describing the method briefly, we use it to generate a complete phase diagram of 
sewing-machine patterns that reproduces all the major features of the diagram determined 
experimentally by Morris et al.'^. We then perform a detailed Fourier analysis of the motion 
of the contact point for each of the patterns simulated, and propose a new classification of 
them based on the spectral content of the motions of the contact point in two orthogonal 
directions. 

For most of the patterns studied, we find that the frequencies present in the spectra of the 
contact point motion are multiples of the steady coiling frequency fic, indicating that the 
dynamics of the sewing machine are closely related to those of steady coiling. Accordingly, 
we set the stage for our study with a brief summary of the latter in the next section. 
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FIG. 2. Steady coiling of a viscous thread, (a) Definition sketch, (b) Coihng frequency as a 
function of height for the parameters of the experiments of Morris et al.'^, calculated numerically 
using a continuation methocP^. The solid curve includes the effect of surface tension (lis = 1.84), 
while the dashed curve is for zero surface tension (lis = 0). The portions of the curves corresponding 
to the gravitational (G), inertio-gravitational (IG), and inertial (I) regimes are indicated. The 
dotted lines show the first three eigenfrequencies of a free viscous pendulum. 

B. Steady coiling 

In steady coiling, the contact point of the thread with the surface moves with a constant 
angular velocity fic along a circle of radius Rc (Figure [2^). In most cases, the thread 
comprises two distinct parts: a long, roughly vertical 'tail' which deforms primarily by 
stretching under gravity, and a helical 'coil' in which the deformation is dominated by 
bending and (to a lesser extent) twisting. Thus the radius of the thread within the coil 
ai is nearly constant. By conservation of mass, the axial speed of the fluid in the coil is 
Ui = Q/Tial 

Steady coiling can occur in several distinct dynamical regimes characterized by different 
balances of the viscous, gravitational, and inertial forces acting on the threacP^^^. These 
regimes appear clearly on plots of the coiling frequency vs. the fall height. For convenience, 
we define a dimensionless fall height H and a dimensionless coiling frequency Cl: 

Figure [2)3 shows ^2 as a function of H for the parameters of the experiments of Morris et al.'^, 
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viz., V = 0.0277 m2s-\ p = lO^kgm-^ 7 = 0.0215Nm-\ d = SmmandpQ = 0.0270gs-^ 
For H < 1.2, coiling occurs in a gravitational (G) regime. Inertia is negligible everywhere in 
the thread, which is governed by a balance between gravity and the viscous forces that resist 
stretching (in the tail) and bending (in the coil). At intermediate heights 1.2 < H < 2.2, 
a complex inertio-gravitational (IG) regime appears, in which the coiling frequency is a 
multivalued function of the fall height. The centrifugal force now becomes important in 
the tail, which behaves as a distributed pendulum with an infinity of eigenmodes whose 
corresponding eigenfrequencies are proportional to the simple pendulum frequency {g/ HY^"^. 
The first three of these frequencies are shown by the dotted lines in Figure [2]d. When one of 
these eigenfrequencies is close to the frequency set by the coil, the tail enters into resonance 
with the latter, giving rise to resonance peaks that appear as rightward-facing 'bumps' in the 
curve Q{H). For large heights H > 2.2, coiling occurs in an inertial (I) regime in which the 
viscous bending force in the coil is balanced by inertia^. The tail in this regime is almost 
perfectly vertical, and is controlled by a balance between gravity, the viscous stretching 
force, and the axial momentum flux. Finally, there is also a viscous (V) regime in which 
both gravity and inertia are negligible everywhere in the thread, but this is only observed 
when both H and d are much smaller than in the experiments of CWL and Morris et al. 

In a typical laboratory experiments on steady coiling, the parameters d and Q and the 
fluid properties p, u and 7 are held fixed while H is varied. Each such experiment is 
completely defined by the values of the three dimensionless groups 

-(^)"'' -(^)"' 

As an example. Hi = 610, 112 = 0.370, and lis = 1.84 for all the experiments of Morris 
et al.^. The (dimensionless) functional dependence of the coiling frequency on the other 
parameters now takes the general form 

Q = Q{H,U,,U2,U,). (3) 

The effect of surface tension is measured by the parameter Surface tension modifies the 
coiling frequency quantitatively but introduces no essentially new dynamics, as can be seen 
by comparing the solid (lis = 1.84) and dashed = 0) curves in Fig. ^p. 

The continuation method used to generate the curves in Fig. [2Jd can be used for steady 
coiling because the flow is stationary in a co-rotating reference frame that moves with the 



contact point. No such reference frame exists for the sewing machine configuration. We 
therefore require a different numerical method, which is described in the next section. 

II. NUMERICAL METHOD 

Our numerical experiments of the sewing machine were set up using the computational 
method of Discrete Viscous Threads (henceforth DVT) originally described in a conference 
papei'^, and presented in detail in a recent papei'^. To the best of our knowledge, DVT 
is the only available numerical method that is fast and robust enough to be applicable to 
the sewing machine geometry while retaining all the relevant modes of deformation, namely 
stretching, twisting and bending. For a thin thread the stretching modulus varies as the 
square of the thread's radius, while the bending and twisting moduli are proportional to 
the fourth power of the radius. As a result the dynamics of thin threads is a nonlinear 
and numerically stiff problem. The DVT method addresses this difficulty by introducing an 
elaborate spatial discretization of the equations based on ideas from differential geometry. 
The method allows simulations to be run for arbitrary mesh sizes, even very coarse ones, 
with optimal stability. This contrasts with conventional discretization schemes which are 
typically stable for sufficiently small mesh sizes only, the maximum mesh size being in 
practice a small fraction of the smallest length scale in the flow, here the small size of the 
coiled region at the bottom. 

Below we briefly present the principles of the DVT method, introduce adaptive mesh 
refinement which provides a tremendous speed-up of the calculations when gravity stretches 
the tail significantly, validate the code against known solutions for steady coiling, explain 
the details of the numerical procedure, and finally present our numerical results. 

A. Smooth case: the equations for thin viscous threads 

The numerical code makes use of the Lagrangian approach and the viscous thread is 
discretized as a polygonal line. A mass is assigned to each vertex, forces are set up on these 
masses, and the motion of each mass is obtained by time integration of the fundamental 
law of dynamics. The discrete forces are designed in such a way that the motion of the 
polygonal line is equivalent to that of a thin viscous thread in the smooth limit. The key 
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element of the numerical method is the expressions for these discrete viscous forces. To 
prepare the derivation, we start by reformulating the smooth case, usually expressed in 
Eulerian variables, in the Lagrangian framework. We introduce a Lagrangian coordinate 
S that marks cross-sections and follows them during motion. This Lagrangian coordinate 
5* is defined as the arc-length in an imaginary reference configuration where the thread is 
a cylindrical tube of constant radius. It plays the same role as the vertex index i in the 
discrete case. 

At a particular time t, the configuration of the thread is defined by its center-line r{S,t) 
and an orthonormal triad {di{S,t),d2{S,t),d3{S,t)). This triad allows one to keep track 
of twisting, because the directions of di and d2 follow the orientation of a cross-section 
as it spins about the center-line. Thin threads deform is such a way that cross-sections 
remain approximately planar and perpendicular to the centerline. As a result, the center- 
line tangent r'{S,t) = || and the unit vector d3{S,t) are aligned. Denoting by i{S,t) the 
axial stretch factor based on the reference configuration, given by the norm of r'{S,t), we 
have: 

r'{S,t)=i{S,t)dsiS,t). (4) 

Since the triad (dj(5', t))i=i_2,3 is orthonormal, its time evolution defines a rigid-body 
rotation for any particular value of S. The associated instantaneous angular velocity uj{S, t) 
is such that 

^^^ = a;(^,t)xd,(5,t) = 1,2,3). (5) 

One can take advantage of the fact that the vectors r' and d^ must remain aligned by 
equation (|4]) to capture the twisting motion of the thread using a single parameter, the 
angular spinning velocity w{S,t) defined by w{S,t) = uj{S,t) ■ d3{S,t). In this view, which 
we call the centerline/spin representatioiJ^, the material velocity lj is a secondary variable 
which can be reconstructed by the equation 

uj{S,t) = ds{S,t) X d3{S,t)+w{S,t)ds{S,t). (6) 

where the time derivative ^ is denoted using a dot. In a viscous thread, the fundamental 
kinematical quantities are the strain rates, defined by 

I dujS^t) l du;iS,t) 
es - J ■ d3(^, t), - - (7) 



where u = r is the center-hne velocity. Here eg captures the strain rate associated with 
axial stretching, while the vector etb captures in a combined manner the strain rates for 
the twisting and bending modes. The strain rate for stretching, ig, is related to the time 
evolution of center-line stretch i by the formula is = d In i/dt. 

For a thin thread, the internal stress is described by the resultant n{S,t) of the viscous 
forces over a particular cross-section, and their moment m{S, t) with respect to the origin 
of the cross-section. These internal force and moment vectors play the same role as the 
tensor cr in 3D continuum mechanics. For the special case of a thin thread geometry, 
Stokes's constitutive law states that stress is proportional to the rate of deformation: for 
the stretching mode, we have 

n(5,t)-d3(5,t) = 3/iAe3, (8) 
and for the bending and twisting modes, 

m(^,t) = (3/i/Pi2 + 2/iJP3) -etb. (9) 

Here fi = pu is the dynamic viscosity, A = ir a"^ the cross-sectional area, / = na^/A is the 
moment of inertia, a is the radius, P3 = ds is the tangent projection operator, and 
P12 = 1 — P3 is the normal projection operator. The radius a{S, t) is a dependent variable 
which is reconstructed by the incompressibility condition 

where = d/2 is the radius in the configuration of reference (note that £ = 1 in the 
reference configuration by definition). The stretching modulus 3fiA was originally derived 
by Trouton!^, and the bending modulus 3/i / can be found for instance irpSl. 

These equations are complemented by the balance of linear and angular momentum, 
known as the Kirchhoff equations for thin rods: 

MM,P(5,,.,..?!fi) (u) 

dm{S,t) dr{S,t) 



dS dS 



xn{S,t) = 0. (12) 



Following an approximation introduced by Kirchhoff himself which is valid for thin threads, 
we have neglected the rotational inertia in the second equation. The vector P{S,t) is the 
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resultant of the external forces per unit reference length dS. The weight of the thread and 
the surface tension are taken into account by setting 

P = Mog- """°'^,f-'^-'" (13) 

Co 

where Aq = it oq^ is the cross-sectional area in the reference configuration, g the accel- 
eration of gravity, and 7 is the surface tension at the fluid- air interface. The last term 



in equation (13) is the net force on the center-line due to surface tension at the fluid-air 
interface, the argument in the derivative being the axial force due to the capillary overpres- 
sure (vra^) (^). Note that there is no need to consider a linear density of applied torque in 



equation (12) for the problem at hand. 

With suitable boundary conditions, the set of partial differential equations (|4-12) con- 
stitutes a well-posed mathematical problem governing the dynamics of a viscous thread. 



B. A variational view: Rayleigh potentials 



The equations of motion (11-12) and the constitutive law (|8j-[9j) can be discretized in 
a natural mannei'^ if they are first re-written in terms of a Rayleigh potential. The 
Rayleigh potential V yields the power dissipated by viscous forces as a function of the 
vertex velocity u(S') = r{S) and spinning velocity w{S) = <-i;(S') • d3(S'). For a thin thread, 
it has three contributions corresponding to the stretching, bending and twisting modes, 
V{u,w) = Vsi^u) + Vt{u,w) + X'b(u,w). Note that the Rayleigh potential V also depends 
on the current configuration r{S,t) but this dependence will be implicit in our notations. 
As an illustration, consider the stretching contribution Pg- It only depends on the vertex 
velocities, and reads 

P,(u) = ^ j 3fiA{e,fidS (14) 

where eg in the integrand is given by equation ([7]), and £dS is the infinitesimal arc-length 
in current configuration. The expressions for the twist and bending contributions can be 
found in referencd^. 

The Rayleigh dissipation potential is useful as it captures the effect of the internal viscous 



stress in a compact mathematical form. Indeed in the equations of motion (11 12) the net 



viscous force n' and the net viscous moment m' + r' x n in the left-hand sides can be shown 
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to be equivalent to a density of applied force 



dV{u, w) 
du{S) 



(u=x(t),ui=ui(t)) 



(15) 



and a density of applied twisting torque 



QAS,t) 



dV{u, w) 
dw{S) 



(u=x(t),u)=ui(t)) 



(16) 



In these equations, the right-hand sides denote functional derivatives, as the dissipation 
potential V takes the functions u{S) and w{S) as arguments. The hat notation expresses 
the fact that the derivatives have to be calculated formally first, and evaluated with using 
the real motion (u = r, w) next. 

C. Discretization 

In the discrete case, the center-line of the thread is represented by a polygonal chain of 
n + 2 particles R(t) = {ro(t), ri(t), ■ ■ ■ , r„+i(t)}. The length t{t) and unit tangent dl{t) of 
an edge i are defined by 



We consider viscous threads having a circular cross-section. As a result there is no need 
to keep track of the absolute orientation of the cross-sections during motion. Twist is taken 
into account through the instantaneous angular velocity of spin of an edge, noted w^{t). This 
is an unknown of the motion, for which we will derive an equation. Representing rotations 
by a single degree of freedom is beneficial for the simulation. The angular velocity vector u}^ 
is a secondary quantity in the simulation, which is reconstructed from the spinning velocity 

by an equation similar to the smooth equation 

The generalized velocity of a viscous thread is obtained by complementing the vertex 
velocities R(t) = {r'o(t), ■ ■ ■ } with the spinning velocities of the edges W(t) = {w^it), ■ ■ ■}. 
The dynamics of the thread is specified by a differential equation involving the position 
R(t), velocities R(t), W(t), as well as the acceleration R(t). Rotational inertia is neglected 
and so W(t) does not enter into the equation. This differential equation is derived next. 

As in the smooth case, internal viscous stress is captured by means of a discrete 
Rayleigh dissipation potential which is the sum of three contributions, '^'(U, W) = Ps(U) + 



r,+i(t)-r,(t)=f(t) dlit). 



(17) 
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'Dt(U, W) + V]^(\J ,W). As an illustration, the stretching contribution is defined in close 



analogy with equation ( 14 ) by 



(18) 



where = 3/i A* is a discrete stretching modulus defined by analogy with equation (|8j), 
and el = ji dy (uj+i — Uj) is a discrete axial strain rate defined by analogy with equation ([t]). 
For an expression of the twist and bending contributions Pt and "Db, we refer the reader 
tcP^. They make use of discrete notions of curvature and twist, based on ideas from discrete 
differential geometry. 

In analogy with the smooth case, we write the equations of motion as 



Pv(t) + P(t) = M-R(t) 
Qv(t) = 0. 



(19) 
(20) 



The first equation is a balance of linear momentum for the vertices and is associated with 
the positional degrees of freedom R(t), while the second equation is a balance of twisting 
torque at each edge, associated with the spinning degrees of freedom W(t). Here, M is the 
diagonal matrix filled with the mass of the vertices, Pv and Qv are the viscous forces and 
twisting torques representing the internal stress in the thread, and P combines the weight 
and the net force on vertices due to surface tension. For a detailed derivation of the discrete 
surface tension forces, se^^. As in equation (12) for the smooth case, we have neglected 



rotational inertia in the right-hand side of equation (20): we have checked that this has 



negligible effect on the simulation when the thread is thin enough. 

Our discretization of the thread is based on the Rayleigh potentials, and the discrete 



viscous forces and moments are defined as in equations (15) and (16) by: 

dV((j,W) 



Pv(t) 
Qv(t) 



dV(U,W) 



(U=R(t),W=W(t)) 



dW 



(21) 
(22) 



(U=R(t),W=W(t)) 

The equations of the present section provide a complete system of equations for the 
dynamics of a discrete viscous thread. For the time discretization, we use a semi-implicit 
Euler scheme, which provides good stability even for quite large time-steps (by semi-implicit, 
we mean that we linearize the equations near the current configuration at every time-step, 
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before applying an implicit scheme). The treatment of boundary conditions is explained in 
section III E[ 

D. Adaptive mesh refinement 

The DVT method uses a Lagrangian grid which is advected by the flow. In sewing 
machine experiments, gravity can typically stretch the centerline by a factor 5 to 10 over the 
course of the descent. In the absence of refinement, this makes the grid very inhomogeneous: 
to capture the small-scale features near the belt one has to use an extremely fine mesh 
size near the nozzle. Overall, a large number of degrees of freedom are required and the 
simulation is inefficient. In addition, important inhomogeneities in edge lengths makes 
the time-stepping problem ill-conditioned and robustness is affected. To overcome these 
difficulties, we have set up adaptive mesh refinement. 

In our implementation of refinement, edges are subdivided whenever their length exceeds 
a nominal length, which is a prescribed function of the distance to the belt. In the upper 
part of the belt, this nominal length is equal to twice the initial segment length, defined 
by the periodic release of new (Lagrangian) vertices from the nozzle. To better resolve the 
coil region, this nominal segment length decreased near the belt according to a prescribed 
exponential profile. This profile was adjusted in such a way that the coil region always 
includes a sufficient number of vertices, typically 10 to 30, with a final edge length typically 
below 0.006 [iP' j g)^^"^ ^ and that the interval between subsequent subdivisions of a given edge 
is always larger than two simulation steps. 

Whenever an edge was marked as needing subdivision, a new vertex was inserted. We 
computed the properties of the new vertex and of the two new edges as follows: the La- 
grangian coordinate of the new vertex is obtained by linear interpolation, the mass stored 
in the former edge is equally split among its children, the position and velocity of the new 
vertex are calculated by an interpolation of order 4, the cross-sectional area A, the spinning 
velocities, and the viscosities of the new edges are obtained by second-order interpolation. 
These interpolation orders were chosen in such a way that the bending and twist forces 
remain continuous upon subdivision. We used the steady coiling geometry to adjust the 
subdivision parameters and to validate the subdivision scheme. 
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E. Emission from the nozzle, capture by the belt 

We found that the implementation of boundary conditions was a key point to successfully 
reproduce the patterns and the phase diagram of the experimental sewing machine. We tried 
simple implementations first, and could reproduce the curves Clc{H) for the frequency of 
steady coiling as well as the simplest stitch patterns, but failed to reproduce entire regions 
of the phase diagram. Further examination revealed the presence of spurious oscillations 
in the calculated acceleration in the steady coihng geometry {V — 0), even though the 
coiling frequency Clc{H) was correctly predicted. We removed these spurious oscillations 
by a more careful account of the boundary conditions both at the nozzle and at the belt, 
as explained below. Suppressing these oscillations appeared to be sufficient to bring the 
numerical predictions in close agreement with the experimental ones. 

New vertices need to be emitted periodically from the nozzle. In a first implementation 
of the clamped boundary conditions there, we considered an infinite string of fiuid particles 
which were moved with the prescribed ejection velocity Q/Aq, until they passed below the 
nozzle and were released. The position of the first vertex clamped inside the nozzle varies 
abruptly as a new vertex is released, and this was the cause of unwanted oscillations. They 
were suppressed by considering a string of two particles in the nozzle, with a fixed position 
relative to the nozzle; the injection velocity is then modelled by steadily increasing the length 
of the second edge, and periodically inserting a new vertex in third position. 

Impact with the belt was first handled by detecting penetration of vertices into the belt 
at the end of every time-step, and constraining their velocity to match the belt's velocity 
at any subsequent time. This also induces large unwanted fiuctuations in the acceleration, 
which can be interpreted by the fact that the vertical momentum resulting from the collision 
is not transferred to the thread until the following time step. The oscillations were removed 
by using a technique known as roll-back. Then, every time-step involves an iteration whose 
aim is to determine the set of vertices undergoing a collision during the time-step: whenever 
an unexpected collision takes place, the step is discarded, time is rolled back, and a new 
time-step is computed, forcing additional vertices to land on the belt at the end of the 
time-step. An additional difficulty in the implementation of roll-back in the context of a 
linearized implicit scheme, is that a only small number of particles can be captured at every 
step. We circumvented this difficulty by using adaptive time refinement. Such an adaptive 
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FIG. 3. Validation of the discrete numerical algorithm for a steadily coiling viscous thread with 
Hi = 610, 112 = 0.37, and II3 = 0. The discrete simulations with the fall height increasing (dashed 
line) and decreasing (solid line) match closely the solution of the steady-state equations obtained 
using an independent continuation methocP^. Arrows denote transient regimes observed in the 
dynamic simulation when the system jumps to a different solution branch after encountering a limit 
point (by contrast the continuation method records a steady but unstable solution, corresponding 
the part of the curve with negative slope). 

time refinement is presumably not needed if a fully (non-linear) implicit scheme is used, 
such as the one presented in Ref.'^. 

F. Validation 

The numerical code was validated by comparing its predictions of the steady coiling 
frequency to the predictions of the continuation method of Ribe^^. The agreement is excellent 
for all fall heights (Figure [s]). The hysteresis of the dynamic simulation in the range 1.1 < 
H < 2.2 is physical, but inaccessible to the continuation method because the latter records 
all steady-state solutions regardless of their stability. 

III. SIMULATION RESULTS 

Our dynamic simulations of the sewing machine patterns, based on the DVT method were 
carried out with non-dimensional quantities. This is achieved by setting the density p, the 
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viscosity /x, and the gravity g to the value 1. This choice imphes that both the length scale 
{iP' I g)^!'^ and the time scale {yjg^^^l^ of the problem introduced in equation are equal 
to 1. The three other physical parameters namely the nozzle diameter, the flow rate and 
the surface tension were chosen to match the values of IIi = 670, 112 = 0.37 and lis = 1.84 
prescribed by Morris et al. experimentgl^^ d = 0.187, Q = 22.9 10~^ and 7 = 1.20 10^1 The 
simulations were initiated from a vertical thread of uniform radius comprising 172 segments 
of equal length, falling from a height H = 0.86 (gravitational regime) onto a belt at rest. To 
avoid dealing with a shock when the thread hits the belt, the simulation was started with 
the bottom end of the thread clamped into the ground. The simulation was run until the 
radius settled to a steady profile as a function of the elevation, and a steady state of coiling 
was established. Next the {H, V) space was sampled by slowly varying or V in turn. 

As an illustration of the capabilities of our numerical technique, we show in Figure |4] 
a simulation of the 'translated coiling' pattern that occurs for relatively low belt speeds. 
Fig. shows a three-dimensional view of the falling thread and the trace it lays down on 
the belt. The simulation time is 0.73 s for one period of this pattern using a 2.6 Ghz Intel 
Core 17 processor and 8 Go of 1067 Mhz DDR3 memory. 

Fig. shows the trajectory of the contact point in the frame of the nozzle (solid line). 
Interestingly, the belt motion breaks the symmetry of steady coiling not only in the longi- 
tudinal x-direction, but also in the transverse ?/-direction. Fig. |4)d also shows the velocity of 
the contact point relative to the moving belt as a function of position along the trajectory 
(thin lines and green arrows). The magnitude of the relative velocity varies significantly over 
a period, in contrast to the meandering pattern for which it is nearly uniform?^. The relative 
speed is maximum at point A where the contact point is moving upstream along the belt, 
and very small at C where the motion is downstream. Finally, fig. ^ shows the amount of 
viscous power dissipated by the various modes, as a function of arc-length along the thread. 
The curves cross each other at a height 2; 0.1 that corresponds to the transition from the 
bending-dominated coil, to the stretching-dominated tail. Thanks to adaptivity, the coil is 
well resolved and the curves for the viscous power dissipation remain smooth there, even 
though they vary rapidly. 

In addition to the dimensionless parameters Hi, 112 and II3 in equation ^ that describe 
the fluid properties and the ejection conditions, the patterns depend on the dimensionless 
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FIG. 4. Simulation of the translated coiling pattern for ITi = 670, 112 = 0.37, Ha = 0, -ff = 0.98 
and V = 0.022 using 182 vertices for 181 segments, (a) Three-dimensional perspective view. A-D 
denote reference points along the thread, (b) Thick gray curve: trajectory of the contact point in 
the frame of the nozzle. The velocity of the contact point relative to the belt is shown at various 
times (thin gray lines) and highlighted when the contact points passes the reference points A-D 
(green arrows) . (c) Viscous power dissipated per unit length dP/ds by the stretching modes (plain 
black and dashed black curves, respectively), and by the twisting and bending modes (thick grey 
curve), as a function of arc-length s, when the contact point passes the apical reference position 



C. All numerical quantities are dimensionless, as explained in the beginning of section III 

fall height H in equation ([l|, and the dimensionless belt speed 

V = V{ug)-^l\ (23) 

Our simulations were carried out by varying H and V for fixed values of Hi = 610 and 
112 = 0.370 corresponding to the experiments of Morris et al.^^. Some of our simulations 
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FIG. 5. Simulated sewing machine patterns (black curves) together with photographs of similar 
patterns observed experimentall}^^. The comparison is qualitative because the parameters used 
in the simulations (Hi = 670, 112 = 0.37, = 0) are not the same as in the experiments oP. 
A: Transient coiling, B: Meanders, C: Alternating loops (figure-of-eight) pattern, D: Disorder, E: 
Double coiling, F: Double meanders, G: W pattern, H: Stretched coiling, I: catenary. (*) Photo 
Courtesy: Chiu- Webster and Lister. 



were done with a surface tension parameter matching the experimental value lis = 1.84; for 
reasons of numerical stability, however, most of the simulations used = 0. 

Figure |5] summarizes all the types of patterns that were encountered when scanning 
the {H, V) plane in the simulations, together with their experimental equivalent^. The 
simulation reproduces nine out of the ten patterns reported by Morris et and observed 
by CWlP. The only missing pattern, the slanted loops, will be discussed later on. To 
ease the understanding, the name 'figure-of-8' is replaced by 'alternating loops' which more 
accurately describes the pattern in Figure [5p. 

In Figure km a phase diagram lists all the numerical patterns encountered in the {H, V) 



19 



space for H < 0.8. The effect of surface tension is included, lis = 1.84 in the simulation. 
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FIG. 6. Phase diagram of the numericahy simulated patterns as a function of dimensionless fah 
height H < 0.80 and dimensionless beh speed V, for Hi = 670, Ha = 0.37 and = 1.84. The 
observations of Morris et al^^ are shown by colored dots for comparison. The typical appearance 
of each pattern (catenary, meanders, alternating loops, translated coiling) is shown in insets. 

For comparison, the patterns observed experimentally by Morris et al.'^ are shown by dots. 
The agreement between the simulations and the experiments is remarkable: all four patterns 
that were observed experimentally in this region of the parameter space, namely translated 
coiling, alternating loops, meanders and catenary, are recovered. In addition the location of 
the boundaries separating the different patterns are reproduced accurately. 

The inclusion of surface tension gives rise to numerical instabilities for heights approxi- 
mately above H > 0.70, which we have not been able to overcome by decreasing the mesh 
size or the time-step. This is why there is no simulation data shown in the lower right 
portion of Figure |6| Since surface tension is not expected to modify qualitatively the dy- 
namics of threads (see Fig. |2]) we investigated the case of larger fall heights without any 
surface tension (Ilg = 0). Five new patterns were observed for larger fall heights, as shown 
in Figure [7} namely double coiling, double meanders, stretched coiling, W-pattern, disor- 
dered pattern, despite the absence of surface tension in the simulations. The new portion 
of the phase diagram, H > 0.8, is very similar to that reported by Morris et al.'^, shown in 
the inset in Fig. [7| In both diagrams, the alternating loops pattern disappears at a critical 
height, beyond which there is a substantial height 'window' containing only simple patterns 
(catenary, translated coiling, and meanders). When the height is increased, three patterns 
having a complex aspect (disordered pattern, stretched coiling and the double meanders) all 
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appear together at nearly the same height. Finally, for some values of the height, disordered 
patterns, shown in grey in the diagram, occur in two separate ranges of the belt speed, with 
stretched coiling in between. 

It is instructive to compare the numerical and experimental phase diagrams with the 
curves Qc{H) of frequency vs. height for steady coiling, calculated with same value of surface 
tension {U^ = for the simulations, II3 = 1.84 for the experiments). These curves are shown 
below each phase diagram in Figure [7j The comparison reveals that some of the more 
complicated patterns (stretched coiling, W-pattern, disordered pattern) appear at heights 
close to that for the onset of the multivalued IG regime in the steady coiling. In the steady 
coiling geometry, it is known that the multivalued regime is caused by the competition of 
several 'viscous pendulum' modes. This suggests that the complex patterns of the sewing 
machine are produced by the non-linear interaction of these pendulum modes. 

Despite their similarities, the numerical (N) and experimental (E) phase diagrams in 
Figure [7] exhibit some systematic differences. In E, double coiling (pink) first appears at the 
same height as the disordered (grey) and stretched coiling (yellow), whereas in N it appears 
at significantly greater heights. Double meanders (purple) have a common boundary with 
the catenary pattern (black) in E, but occur only for significantly lower belt speeds in N. 
In N, the catenary state can transition to disorder (grey) over a range of heights, unlike 
in E. The range of belt speeds for double coiling is significantly wider in N than in E. 
Finally, in N the W-pattern is observed sporadically and for greater heights than in the 
diagram in Fig. [7} Some of these differences are due to the absence of surface tension in 
the simulations, and to the fact that collisions of the viscous thread with its trace are not 
prevented. Another explanation for the discrepancies may be the fact that Morris et al.'^ 
performed their pattern recognition visually, whereas we used a more quantitative automatic 
procedure based on Fourier decomposition. This is described in the following section, where 
we review each of the individual patterns in detail and propose a systematic classification 
of them. 

IV. ANALYSIS OF THE PATTERNS 

To illustrate our pattern analysis, we fix the belt speed V = 0.02 and vary the fall height, 
thereby following the horizontal line AB in the phase diagram of Fig. [7j In order of increasing 
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FIG. 7. Phase diagram of sewing machine patterns determined from numerical simulations with 
Hi = 670, 112 = 0.37 and no surface tension (lis = 0). The patterns simulated include translated 
coiling (red), meanders (blue), alternating loops (green), double coiling (pink), stretched coiling 
(yellow) , double meanders (purple) , and disordered patterns (grey) . The catenary pattern in black 
would extend all the way up in the diagram, the border between black stripes and the white region 
at the top being only the result of the simulation being stopped for large enough values of V . The 
W-pattern (yellow dots circled in black in inset) and slanted loops (blue dots circled in black in 



inset) are discussed in section IV B The horizontal path AB is used to construct Figure ^ The 



coiling frequency Vt for steady coiling is shown as a function of height H below the phase diagram. 
Inset: Phase diagram determined experimentalljff^ for Hi = 670, 112 = 0.37, and lis = 1.84 (top) 
together with the corresponding curve Q.{H) (bottom). 
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height, the patterns seen along this hne are meanders, alternating loops, translated coiling, 
and double coiling. We track the spectral content of these patterns continuously as they 
change smoothly or bifurcate. To do so, we focus on the trajectory of the contact point of 
the thread with the belt. Let x{t) and y(t) be its longitudinal and transverse coordinates 
in the laboratory (nozzle) reference frame, and define x(t) = x{t) + iy{t). Let X{t,t*) and 
Y{t,t*) be the coordinates (also relative to the nozzle) at time t of a material point that 
was laid down on the belt at time t* < t, and let X{t, t*) = X{t, t*) + i Y{t, t*) be a generic 
point in the trace. The advection by the belt is expressed by: 

X{t, t*) = x(t*) + (t - t*) V. (24) 

This equation means that the pattern X{t, t*) seen on the belt is obtained by unfolding the 
motion of the contact point motion x{t), as illustrated in Fig. [9^. 

The numerical traversal of the line AB in Figure [7] required about 78500 time steps of 
size 6t = 0.1{iy/g'^y^^. We performed a Fourier analysis of the motion x{t) over a sliding 
window of 2000 time steps (the FFT spectrum was computed every 500 steps). The spectra 
obtained in this way typically comprise several well-defined peaks whose frequencies can be 
determined precisely (Fig. |9}d). 

Let u}[^\ UJ2^\ etc. be the peak frequencies of the motion in the x-direction, and u}^\ u}^\ 
etc. be those for the motion in the ^/-direction. Because the fall height H is slowly changing 
with time during the simulation, each observed frequency uj^^^ or can be plotted on 
a diagram to provide a 'portrait' of the evolving frequency content of the contact point 
motion, as a function of the fall height H. The result is shown in Figure [sj The principal 
observed frequencies Un^ and Un '' are indicated in red and blue, respectively. Also shown 
for reference is the steady coiling frequency Qc{H) for the same fluid properties and ejection 
parameters (solid line), together with several multiples (1/3, 2/3, 1/2, 4/3, 5/3, 2) of that 
frequency (dashed lines). 

The first pattern (0.62 < H < 0.8) is meandering, which is characterized by two fre- 
quencies with a ratio ui[^^ /uj^^ = 2. At the lowest height H = 0.62 where meandering first 

(x) 

appears, the meandering frequency uj\ is very close to the steady coiling frequency fic for 
the same height, as expected on theoretical ground^^^*^. The next pattern (0.8 < H < 0.9) is 
the alternating loops, which has a rich spectrum involving five multiples of f2c/3. Translated 
coiling appears next (0.9 < H < 1.35), and is characterized by a single frequency 
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FIG. 8. Frequency content of patterns encountered along the line AB through the phase diagram of 
Figure [Tj Frequencies of transverse (y-direction) and longitudinal (x-direction) motion are shown 
as blue and red, respectively. Note that these frequency are the same in the upper region of the 
diagram where the two colors overlap. In the dashed insets, the patterns are identified using the 
same color codes as in Figure [7j Grey bands indicate ranges of heights for which the patterns are 
transient. Also shown is the frequency of steady coiling as a function of height (brown solid 
line), together with several multiples of that frequency (brown dotted lines). 

very close to the steady coiling frequency. Finally, double coiling (1.35 < H < 1.5) has two 
frequencies u^^^ = u^i ^ ^ Qc and uj2^^ = ^ 

Figure |8] shows that the stitch patterns are combinations of motions in two orthogonal 
directions with frequencies closely related to the steady coiling frequency Cl^. Accordingly, 
we now change our perspective and classify the patterns based on their frequency content 
rather than on the shape they lay down on the belt. This frequency analysis is used to set 
up an efficient tool for identifying the patterns and assembling the numerical phase diagram 
automatically. In addition it leads to a simple kinematic model that provides a unified 
description of the different patterns. 
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Numerical Simulations 




FIG. 9. Three-stage analysis of the patterns, illustrated for the case of meandering, (a) From the 
trace of the thread on the belt (al), we extract the trajectory of the contact point in the frame of 
the nozzle (a2). (bl) Fourier spectra of the longitudinal (red) and transverse (blue) components 
are obtained by Fast Fourier Transform, and compared to the frequencies Q of steady coiling for 
the same height: the blue and red dots in (b2) are the main frequencies in the simulation, while 
the intersection of the horizontal line with the thick black curve is the frequency of steady coiling. 

In the example shown here the belt speed is close to the critical value for the onset of meandering 

(v) 

and the dominant transverse frequency uj{ is close to the steady coiling frequency ilc for the same 



height, (c) The two main frequencies extracted by FFT are injected into the kinematic model (26) 
to generate a synthetic motion for the contact point in the nozzle frame (cl), and a synthetic stitch 
pattern (c2). In the example shown, the motion involves only one transverse frequency uj^^ and 
one longitudinal frequency u)^^ = 2uj^\ 
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A. Spectral signature of a pattern 



We illustrate our method using the example of meandering which in most cases is the first 
pattern to bifurcate from the catenary state as the belt speed decreases. The red and blue 
lines in figure [9]d show typical spectra of the motion of the contact point in the longitudinal 
X and transverse y-directions, respectively. The amplitude of the transverse motion is much 
greater than that of the longitudinal motion, and the frequency of the latter is exactly twice 
that of the former, uj[^^ = 2uj^\ This suggests that a meander pattern can be synthesized 
by retaining only the two main frequencies, viz. 

x{t) = ai cos(2a;i^^t) + i (3i cos{J^^ + 7r/4) (25) 



where ai and /3i are the amplitudes of the longitudinal and transverse motions, respectively. 
Here the phase difference 7r/4 is required to reproduce the symmetry of the pattern. A similar 
two- frequency model was used by Morris et al.'^ to analyze weakly non-linear meanders. 
Figure shows the contact line trajectory in the frame of the nozzle implied by (25) 



with ai/ Pi = 0.2, and Figure |9p2 shows the corresponding meander pattern obtained by 
advecting the motion (25) in the x-direction with V = 1.4(z/(yf)^/^ and = 1. Based on 



Figure |9| we define as a 'meander' any pattern whose longitudinal motion, compared to the 
transverse motion, has twice the frequency, a much smaller amplitude, and a phase shift of 
7r/4. 

Generalizing the above example, we show that all the sewing machine patterns can be 
represented by a superposition of a few harmonic motions of the form: 

x{t) + iy{t) = J2 "i cos{ujp + 0J) + ^ XI cos(wJt + (26) 
j=i i=i 

Two terms in each direction are sufficient to reproduce the main features of the patterns. 



Nx < 2 and Ny < 2. In Eqn. (26), aj and /3j are the amplitudes of the components of 
the motion with frequencies Uj^^ and u}^^\ and (pj and 0j are the phases relative to the 
highest-frequency mode. We now show how each of the sewing machine patterns can be 



characterized in terms of the parameters that appear in Eqn. (26). 
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B. Identification of the patterns 



The identity of each pattern is determined not by the absolute values of the parameters 



in Eqn. (26), but rather by the dimensionless groups that can be formed from them, namely 
frequency ratios, amplitude ratios, and the relative phases 0j. In the following we identify 
the characteristic values of these groups for each of the patterns in turn. For ease of refer- 
ence, the results are summarized in Table [1} 



1. Translated coiling. This pattern occurs for 0.5 ^ -ff ^ 1.35 and low belt speeds 
(Figure [T]). Figure 10 A. shows a simulation of this pattern (upper left) and the correspond- 
ing Fourier spectra of the longitudinal and transverse motions of the contact point (upper 
right), and the steady coiling frequency Clc- The longitudinal and transverse motions have 
similar amplitudes and are characterized by a single dominant frequency that 
is very close to the steady coiling frequency Clc{H) for the same fall height H. The peak 
frequency deviates from its original value as the belt speed increases. The amplitudes in 
both directions remain equal and an almost circular shape is created. The panels at lower 
left show the reconstructed motion in the frame of the nozzle (right) and on the belt (left), 
calculated using (26) with = Ny = 1, ai = Pi = 1 and uj[^^ = oo^^ = 1. Note that the 
experimental pattern shifts upwards or downwards as the belt speed is increased; this shift 
does not affect the spectral context, but could be taken into account by including a purely 



imaginary constant in Eqn. (26). 



2. Meanders. On Figure [7| this pattern is seen for 0.6 ^ ^ 1.3 and a range of 
intermediate belt speeds. The typical Fourier spectra of the meandering pattern were pre- 
viously shown in Figure |9| The pattern is a superposition of a single longitudinal, and a 
single transverse harmonic motions, with frequency ratio uj[^^ /uj^^ = 2, amplitude ratio 
1; and a relative phase = 7r/4 (with the convention (pf = 0). Near the cate- 
nary/meander boundary uj^^ ^ 0,c', farther from the boundary, the meandering frequency 
departs significantly from Clc. The regular symmetrical meanders correspond to a phase 
difference of 7r/4 between the two directions (fig. |9| However, the pattern may deform into 
a bean-like shape in certain cases. This situation was reproduced kinematically by reducing 
the phase difference to a value close to 7r/6. 
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FIG. 10. Kinematic analysis of individual sewing machine patterns. (A) coiling; (B) alternating 
loops; (C) double coiling. For each pattern, the upper left, right, and lower left panels correspond 
to parts a), b) and c), respectively, of Figure |9j The kinematical reconstruction is a proof of 
concept, there is no attempt to match the wavelength of the simulations. 
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TABLE I. Parameters of the kinematic model in Eqn. (26) used to construct synthetic patterns. 
Because the patterns are defined by the relative values of the frequencies and uj\, the frequency 
1 is assigned by convention to the peak of largest amplitude. The other frequencies are given by 
ratios of small integer numbers. A star indicates a frequency that is locked to the steady coiling 
frequency Ctc- Likewise, the amplitudes qi, 02, /3i and (32 are given relative to each other, and 
correspond to typical values. V indicates the speed of the belt used to unfold the synthetic patterns. 
A dash indicates that the parameter is not relevant for the pattern in question. The disordered 
pattern is not reconstructed using the kinematical model as it involves more harmonics. 



3. Alternating loops. This pattern was called 'figure-of-eight' by CWLP and Morris et 
al.'^, but we prefer to call it 'alternating loops'. The domain of this pattern is an elongated 
'bubble' sandwiched between translated coiling and meandering at relatively low fall heights 
H < I (Figure [7|. This pattern displays a remarkably rich frequency spectrum with five 
principal peaks (Figure [lO^-c). In contrast to meandering, the motion with the largest 
amplitude is longitudinal, with a frequency u[^^ that locks onto the frequency 2Qc/3 (see 
also Figure [s]). The next largest peaks correspond to a transverse motion with frequencies 
= Clc and = Clc/3, both amplitudes being very close. The harmonics 4fic/3 in the 
longitudinal direction, and 5Clc/3 in the transverse direction are also visible. 
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The frequencies of all five peaks can be written compactly as 

3 



uj^y) = {2p-l)^ (p= 1,2,3), 



Even though the spectra in Fig. lOB-b shows five peaks, the kinematic model in Eqn. (26) 
generates an almost identical pattern if one retains only the three main contributions u\ , 
uj'f\ and with an amplitude ratio Pi/ai ~ .5, /3i/a2 ~ -5, and phases 0f = 7r/2 and 
(f)\ = (f)2 = (Figure [To^-c). We used these characteristics of the three main frequencies as 
a criterion for automatic detection of alternating loops. 



4., 5. Double coiling and double meanders. The Fourier spectra for these patterns are 
shown in Figures [TO| G and [TI}\, respectively. Both the longitudinal and transverse compo- 
nents have peaks at two frequencies = = Qc/'2 and = u^2 ^ = Qc- The origin 
of these frequencies is clear from the uppermost part of Figure [8} which shows them as 
functions of fall height for double coiling. The range of fall heights in question is within the 
inertio-gravitational regime of steady coiling, for which the curve Clc{H) is multivalued in 
specific height ranges (Figure [2]). The portion 1.2 < H < 1.5 of that curve has two stable 
branches corresponding to different 'pendulum' modes of the tail: a lower branch (labelled 
LB in Figure [sl with ilc ~ 1.15, and an upper branch (UB) with Clc ~ 2.1-2.2. Figure [s] 
shows that the higher double coiling frequency Wg^^ = 00^^ stays locked to the upper branch 
of Clc{H), which is the only stable one when H > 1.37. The lower frequency, by contrast, 
follows a 'phantom' branch with frequency Clc/2 that is very nearly a direct continuation 
of the lower branch of Clc{H) to greater fall heights. This behavior is possible because the 
ratio of the frequencies of the upper and lower branches happens to be quite close to 2.0. 

Although double coiling and double meanders have the same frequency content, they 
are distinguished by the relative amplitudes and phases of the transverse and longitudinal 
motions. For double coiling, the amplitudes of the two motions are the same at both fre- 
quencies («! = /?!, a2 = (32), and the relative phases are ^'j^ = 02 = and (f)\ = (p^ = 0. 
For double meandering, by contrast, the transverse motion is dominated by the frequency 
= Qc while the longitudinal motion is dominated by u^'^ = Qc/2. The relative phases 
are (pf = n/A while = 0. 
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FIG. 11. Same as Figure 10, but for (A) double meanders, (B) stretched coiling, and (C) the 
'W-pattern'. 6 A denotes the difference of amplitude. 



6., 7. Stretched coiling and the W-pattern. These patterns occur predominantly in the 
range of heights corresponding to inertio-gravitational coihng (right-hand side of Figure [?]). 
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Their typical Fourier spectra are shown in Figure [Tip and C, respectively. Like double coil- 
ing and double meanders, their characteristic signature is the {CLc/2,CLc) frequency couple. 
But whereas double coiling and double meanders are dominated by transverse motion at the 
frequency Vtc, stretched coiling and the W-pattern are dominated by longitudinal motion 
at the frequency VLc/2. The difference between stretched coiling and the W-pattern is only 
due to different relative amplitudes of the transverse and longitudinal motions denoted 5a 
in Figure [ll^-b and C-b. The a;— direction is dominant in both cases but the peak at tic 
in the ?/-direction is much closer for the W-pattern than for the stretched coiling. Theses 
differences cause the change of invagination between the two patterns (Figure [TTj3-a and 
C-a). The phase difference between the longitudinal and transverse motions is 7r/2 in both 
cases. 



8. Disorder. Disordered patterns appear in several parts of the phase diagram (grey in 
Figure [7]), primarily at heights within the inert io-gravit at ional coiling regime. The typical 
Fourier spectra of these patterns is very rich, with more than four peaks in both the longi- 
tudinal and transverse directions with comparable and strongly time-dependent amplitudes 



(Fig. 12 -i and corresponding FFT). Such patterns are not transient between two steady 
patterns, as the aperiodic behavior persists indefinitely in time. 



9. Catenary. The catenary is obtained when the point of contact is at rest in the nozzle 
frame, which happens in the upper region of the phase diagram. The FFT spectrum is then 
empty. 



In addition to the patterns discussed above, slanted loop^ were also reported by Morris 
et al.'^ in a very narrow region of their phase diagram. Slanted loops is a pattern wherein 
a buckle is periodically laid down on the belt, and subsequently closes up into a loop when 
the thread touches itself and coalesces. Our numerical simulations do not account for self- 
contact of the thread, nor for surface tension-mediated coalescence. This probably explains 



why we observed slanted loops transiently only, as shown in Fig. 12 -ii. A similar argument 
may also explain why the W-pattern is observed for significantly larger fall heights in the 
simulations than in the experiments: considering self-contact of the thread would certainly 
favor its existence over the stretched coiling pattern in the simulation. 
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i) Disorder 




t + At 



ii) Slanted loops 



t + At 




FIG. 12. Unsteady patterns: (i) disordered pattern and (ii) slanted loops. Both patterns are shown 
at two different times of the simulations, and compared in the insets to experiment^. The lower 
part of the figure illustrates the FFT of the disordered patterns. 



We observe that CWlP reported yet another pattern, side kicks. Side kicks consist of small 
heaps of fluid regularly spaced along an otherwise perfectly straight trace. We suggest that 
this pattern is a limiting case of stretched coiling in which the amplitude of the transverse 
motion becomes very small relative to that of the longitudinal motion. Side kicks have not 
been reported in the experimental phase diagram of Morris and al., the one we attempt to 
reproduce in the simulations; consistently, we have observed this pattern in our simulation 
only transiently. 



V. DISCUSSION 

The simulations of the fluid-mechanical sewing machine presented here were performed 
using a new numerical algorithm. Discrete Viscous Threads (DVT). The essential idea of 
DVT is to start from a complete geometrical and kinematical description of the thread in the 
discrete setting, and push the discrete approach as far as possible; in particular a discrete 
representation of viscous stress is built based on a variational view (Rayleigh potentials). 
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This approach leads to a code that is robust even for quite large mesh size. This numerical 
method is used to simulate complex unsteady behavior of a viscous thread and offers good 
efficiency. As an example, the traversal (78500 time steps) of the line AB in Figure [T] required 
987 s on a 2.87 GHz Intel Core processor. The accuracy of the method is demonstrated by 
its ability to reproduce the curve of steady coiling frequency vs. height, as predicted by 
an independent continuation method (Figure [2]), and by the close agreement between the 
calculated and experimentally determined phase diagrams for the sewing machine patterns 
(Figure |6]) . 

At each time step in an unsteady DVT simulation, any desired kinematical or dynamical 
variable can be calculated as a function of arc-length, i. e. at every vertex along the thread's 
centerline. Examination of these functions provides insights into the thread's dynamics. 
We saw an example in Figure which showed the local rates of viscous dissipation of 
energy due to deformation by stretching (or compression), bending, and twisting. The figure 
reveals that the thread is divided into two distinct parts: a 'tail' in which the dissipation 
is dominated by stretching, and a 'coil' in which it is dominated by bending and twisting, 
but with a significant contribution from compression. The differential equations describing 
bending are of higher order than those describing stretching, so Figure ^ implies that the 
coil is an 'inner' solution or boundary layer whose presence is required by the need to satisfy 
all the relevant boundary conditions at the thread's contact point with the belt. While the 
boundary-layer character of the coil region has long been recognized for steady coilin^^^*^, 
our simulations open up the possibility of studying the associated non-steady dynamics. 

Another benefit of the simulations is to allow exploration of regions of parameters space 
that are inaccessible in the real world but provide new insights into the dynamics of the 
thread. In Figure the rates of viscous dissipation for the bending and twisting modes 
were added together. A real thread is made up of an incompressible fluid, and the ratio of the 
bending to the twisting modulus is always 3/2 (both moduli are in particular proportional 
to the fourth power of the thread's radius). In the simulation, it is possible to investigate 
the relative importance of bending and twisting in the thread's dynamics by setting the 
twisting modulus to zero, taking 2/i/ = in equation ([9]), while keeping the bending mod- 
ulus unchanged. We performed additional DVT simulations for a perfectly twist-compliant 
sewing machine. The resulting phase diagram shows only minor differences with Figure [7], 
showing that twist plays a negligible role, relative to bending, in the selection of the stitch 
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patterns. 

The experiments and simulations of the viscous sewing machine reveal in the first instance 
a great diversity of patterns, whose relations to one another are not evident. Our goal was 
to characterize the patterns in a more unified way. This is possible by going beyond a visual 
identification and computing for each pattern the Fourier spectra of the longitudinal and 
transverse components of the motion of the thread's contact point with the belt. We showed 
that each pattern has a distinct spectral signature consisting of isolated peaks at a small 
number of well-defined frequencies. The patterns differ from each other in the values of those 
frequencies, in their relative amplitudes, and in their distributions among the longitudinal 
and transverse modes. 

A closer look shows that the frequencies in the spectra are closely related to the frequency 
Clc of steady coiling of a thread falling on a motionless {V = 0) surface from the same height. 
The precise nature of the relationship depends on the pattern considered. For meanders, the 
frequency of the transverse motion ^ Qc at onset, but then deviates significantly from Qc as 
the belt speed is decreased beyond the critical value (Figure [s]). In all the other patterns, 
however, the frequencies are locked to Qc in some way. In stretched coiling, the dominant 
frequency of both the longitudinal and transverse motions is Cl^. Still more complicated are 
double meanders, double coiling, stretched coiling, and the W-pattern, for which the dom- 
inant frequencies are Clc and Clc/2. The presence of these frequencies reflects the nonlinear 
interaction of the two lowest modes of inertio-gravitational coiling, whose unforced frequen- 
cies differ by approximately a factor of 2 (Figure [s]). Among periodic patterns, the richest 
spectral content is achieved by the alternating loops pattern, for which the five dominant 
frequencies are multiples of (Figure [8|. 

We proposed a simple kinematic model whereby each pattern is reconstructed by a su- 
perposition of a few frequencies, with appropriate amplitudes and relative phases in the 
longitudinal and transverse directions. This makes it possible to set up automated recog- 
nition of the patterns, and leads a classification of the patterns within a unified descriptive 
framework. The next step is to elucidate the physical mechanisms responsible for these 
simple spectral signatures. In future work we plan to investigate the similarities between 
the sewing machine and low-dimensional oscillator models with non-linear forcing, using the 
direct DVT simulations as a starting point to look into these analogies. 

We are very grateful to Miklos Bergou and Eitan Grinspun as the collaborative devel- 
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opment of the DVT method has made the present work possible. Prehminary numerical 
results concerning the viscous sewing machine have been obtained in collaboration with 
themfHEB would like to thank E. Grinspun for suggestions on the manuscript. During 
the preparation of this manuscript, we learnt of an independent but related experimental 
work by R. L. Welch, B. Szeto, and S. W. Morri^. 
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